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ABSTRACT 

This  report  discusses  the  results  of  an  extensive  long-range  reversed  refraction 
profile  which  traverses  the  Michigan  Basin  and  the  Appalachian  Mountains.  Particu¬ 
lar  emphasis  is  placed  on  the  attenuation  of  first  compressional  wave  arrivals  and  on 
crustal-upper  mantle  ntructure.  Also  included  are  results  of  a  survey  on  lake  bottom 
seismic  background  noise  in  Lake  Superior.  Although  large  signal  levels  were  ob¬ 
tained  in  the  lake  bottom,  long  periods  of  high  background  noise  indicated  that  land 
based  seismographs  are  superior.  Theoretical  studies  in  elastic  wave  scattering 
were  made  to  evaluate  the  effect  of  nonhomogeneities  on  the  propagation  of  elastic 
waves.  Theoretical  development  of  mode  filtering,  (previously  demonstrated  empiri¬ 
cally),  is  nearly  completed.  A  LASA  digital  data  convertor  was  designed,  constructed, 
and  put  into  operation.  Two-  and  three-dimensional  seismic  ray- tracing  techniques 
using  digital  modeling  were  developed.  These  techniques  enable  the  seismologist  to 
determine  the  effects  of  lateral  inhomogeneities  and  irregularities  in  velocity  inter¬ 
face  on  travel  time  and  attenuation  of  seismic  waves. 
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COLLECTION  AND  ANALYSIS  OF  SEISMIC- WAVE  PROPAGATION  DATA 

Annual  Report 

1  June  1967  Through  31  May  1968 
SUMMARY 

An  extensive  long-range  reversed  refraction  profile  was  obtained  in  the  eastern  United 
States.  This  profile  traversed  the  Michigan  Basin  and  the  Appalachian  Mountains.  Close  sta¬ 
tion  control  (25  km  or  less)  over  most  of  the  profile  permitted  detailed  studies  of  travel  time 
and  attenuation  anomalies.  Crustal  and  upper-mantle  structure  was  determined  along  the  pro¬ 
file.  A  mountain-root  system  similar  to  the  Rocky  Mountains  was  indicated  under  the  Appala¬ 
chian  Mountains.  A  two- layered  crust  with  a  high-velocity  (7.7  km/sec)  lower  layer  was  found 
in  the  Lake  Superior  area,  and  a  single-layered  crust  was  determined  in  the  Atlantic  Coast  area. 
High  upper-mantle  Q  values  (Q^  a  1900)  for  first  compressional  wave  arrivals  were  disclosed 
in  the  Lake  Superior  area  and  across  the  Michigan  Basin;  considerably  lower  values  (Qq  s  960) 
were  indicated  farther  east.  Lower  crustal  Q's  (Q^  s  160  and  200)  were  found  for  the  two-lay¬ 
ered  crust  in  the  Lake  Superior  area;  a  considerably  higher  value  (Qa  a  870)  was  determined 
for  the  single-layered  crust  along  the  East  Coast. 

A  new  method  of  displaying  seismic  attenuation  and  spectral  data  was  developed;  it  permits 
the  presentation  of  the  data  in  three  dimensions  (distance,  frequency,  and  ground- particle  veloc¬ 
ity  amplitude).  Meaningful  and  consistent  anomalies  in  the  spectrum  of  seismic  waves  were 
found  that  could  be  correlated  with  variations  in  the  source  and  in  propagation  effects. 

Theoretical  studies  in  elastic-wave  scattering  were  made  to  evaluate  the  effect  of  non¬ 
homogeneities  on  the  propagation  of  elastic  waves.  A  general  solution  was  obtained  to  determine 
the  scattering  of  plane  waves  in  an  elastic  medium  by  elliptical  obstacles.  The  degenerate  case 
of  a  slit  (fault)  was  also  computed.  Current  studies  are  directed  toward  evaluating  the  scatter¬ 
ing  effects  caused  by  imperfections  or  irregularities  at  a  velocity  discontinuity  on  the  reflected 
and  transmitted  seismic  waves. 

A  series  of  seismic  recordings  was  obtained  on  the  bottom  of  Lake  Superior  to  determine 
the  nature  of  the  background  noise  and  to  check  the  reciprocity  in  signal  amplitudes  of  tele- 
seismic  sources.  Previous  research  had  disclosed  that  shots  fired  in  Lake  Superior  were  very 
efficient  in  the  coupling  of  seismic  energy.  If  teleseismic  sources  recorded  on  the  bottom  of 
Lake  Superior  had  large  amplitudes  and  the  background  noise  were  low,  consideration  might  be 
given  to  establishing  a  permanent  station  or  array  for  monitoring  small  distant  events.  The 
results  of  this  investigation  disclosed  that  larger  amplitude  signals  were  obtained  on  the  lake 
bottom  but  that  the  signal-to-noise  ratios  were  no  better  and  in  many  cases  worse  than  those 
obtained  on  land  stations.  Long  periods  of  high  background  noise  were  encountered  on  the  lake- 

bottom  seismograph  and  could  be  correlated  with  local  meteorological  conditions  and  with  water- 
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wave  action. 


Analysis  of  travel-time  and  amplitude  data  from  the  SHOAL  event  at  distances  of  200  to  1200 
km  has  revealed  significant  azimuthal  variation  in  velocity  and  attenuation.  Pn  travel- time  re¬ 
siduals  east  of  the  Sierras  are  negatively  correlated  with  teleseismic  P  residuals  and  can  be 
explained  in  terms  of  rather  localized  velocity  variations  in  the  uppermost  mantle.  Average 
apparent  Q  values  were  276  for  Pn  and  225  for  Pg.  Higher  Q's  were  observed  to  the  north  for 
Pn  and  to  the  southeast  for  P^.  These  results  indicate  that  explosions  detonated  in  areas  of 
complex  crust  and  upp'r-mantle  structure  may  be  expected  to  radiate  seismic  energy  in  a  high¬ 
ly  asymmetrical  manner. 

The  mode- filtering  techniques  developed  at  the  Geophysics  Laboratory  by  optical  processing 
methods  have  moved  from  an  empirical  basis  to  a  nearly  complete  theoretical  basis. 

Given  the  sampling  geometry  described  in  the  previous  annual  report  under  this  contract  [1], 
it  can  be  shown  that  filters  which  are  merely  opaque  and  transparent  quadrants  symmetrically 
placed  with  respect  to  the  origin  yield  an  output  which  is  the  sum  of  the  original  seismic  traces 
and  those  whose  phase  spectrum  has  been  shifted  by  tt/2.  The  use  of  banded  filters  of  the  same 
symmetry  as  above  provides  (no  proof  has  been  given  yet)  amplitude  equalization  before  (or 
along  with)  the  phase  shift  to  further  attenuate  Rayleigh  waves.  Of  course,  the  phase  shifting 
must  not  be  applied  to  signal  components;  they  must  be  identified  in  the  Fourier  plane  and  the 
liters  applied  only  to  noise  patterns  in  that  plane.  To  dace  this  technique  has  been  directed  to 
ward  extracting  P-wave  signals  out  of  Rayleigh  wave  noise  where  signal-to-noise  ratios  are 
frequently  less  than  unity.  Plans  are  being  made  to  extend  the  technique  to  include  S-wave  sig¬ 
nal  enhancement. 

A  digital  data  convertor  for  LASA  data  was  designed,  constructed,  and  is  operational.  It 
permits  viewing  up  to  50  data  channels  directly  from  a  LASA  digital  tape  without  the  use  of  a 
digital  computer.  The  convertor  is  a  flexible  device  with  either  an  analog  or  a  digital  output. 
Wiggly  line  seismograms  or  variable- density  film  records  can  be  readily  obtained.  An  investi¬ 
gation  of  on-ll-.e  optical  techniques  (velocity  filtering,  beam  forming,  cross  correlation,  etc.) 
for  the  detection  of  significant  events  is  now  under  consideration.  While  the  system  would  not  be 
as  versatile  and  as  accurate  as  the  elaborate  digital  facilities  at  SAAC,  the  technique  may  provide 
a  low-cost  means  of  viewing  LASA  data  in  a  routine  manner  and  aid  in  the  selection  of  events 
for  more  elaborate  processing  by  a  digital  computer. 

A  three-dimensional  seismic-ray-tracing  technique  and  travel-time  computation  using  digi¬ 
tal  modeling  is  under  development;  this  will  enable  the  seismologist  to  determine  the  effects  of 
lateral  inhomogeneities  and  irregularities  in  velocity  interfaces  on  the  travel  times  and  ampli¬ 
tudes  of  seismic  waves  emerging  at  the  surface.  The  geological  model  in  any  or  all  portions  of 
the  wave  path  from  the  source  to  receiver  can  be  varied,  including  the  underlying  structure  of  an 
array.  It  is  anticipated  that  simulation  of  the  type  of  anomalies  observed  at  LASA  will  be  possi¬ 
ble.  The  preliminary  results  of  this  investigation  have  been  quite  encouraging.  Calcomp  plots 
of  the  digital  computer  output  are  shown  for  two-and  three-dimensional  models.  The  output  will 


also  be  displayed  on  a  CRT  to  allow  the  interpreter  to  rapidly  vary  the  geological  structure. 
Also,  iterative  perturbation  of  the  velocity  structure,  planned  for  the  computer,  will  alter  the 
velocity  structure  until  travel-time  data  is  duplicated. 
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SEISMIC- WAVE  PROPAGATION—  TRAVEL  TIME  AND  ATTENUATION 

D.  E.  Willis 

1.1.  INTRODUCTION 

During  the  past  five  years,  The  University  of  Michigan  has  participated  in  a  number  of 
cooperative  investigations  of  the  earth's  crust.  A  large  portion  of  this  work  involved  a  long- 
range  reversed  seismic  refraction  profile  extending  along  a  line  from  International  Falls,  Min¬ 
nesota,  to  Jacksonville,  North  Carolina,  a  distance  of  approximately  2000  km.  The  shot  location 
and  recording  stations  are  shown  in  figure  1.  The  major  shot  programs  include: 

1962  North  Carolina  Experiment 

1963  Lake  Superior  Experiment 

1964  July  Lake  Superior  Experiment 

1964  October  USGS  Shots 

1965  ECOOE 

1966  EARLY  RISE 

1965  CHASE  m,  700  tons 

1965  CHASE  IV,  310  tons 

1966  CHASE  VH,  400  tons 

The  section  of  the  map  shown  here  is  from  a  Calcomp  plot  made  through  use  of  an  IBM  7090 
computer.  All  of  the  geographical  features,  the  shot  points,  and  recording  stations  were  machine 
plotted. 

In  this  report  we  present  some  of  our  data  on  reduced  travel  times,  record  sections,  and 
attenuation,  correlating  as  much  of  them  as  possible  with  the  geological  environment. 

1.2.  REDUCED  TRAVEL  TIMES 

In  figure  2  the  reduced  travel  times  are  shown  for  the  1-  and  5- ton  shots  recorded  on  the 
ECOOE  northern  profile.  An  upper- mantle  velocity  of  8.2  km/sec  was  used  for  these  computa¬ 
tions.  There  appears  to  be  considerable  scatter  on  this  graph.  However,  if  each  station  were 
separated,  it  could  be  seen  that  the  stations  located  east  of  the  Appalachians  group  together, 
giving  an  intercept  time  of  approximately  +7  1/2  sec,  and  those  over  the  crest  and  just  west  of 
the  Appalachians  group  together  with  an  intercept  time  of  approximately  9  to  9.5  sec.  This 
could  be  caused  by  a  mountain- root  system  similar  to  that  of  the  Rocky  Mountains. 

Figure  3  shows  the  reduced  travel-time  data  for  the  1-  and  10-ton  shots  of  the  ECOOE 
southern  profile  and  the  CHASE  III  and  VII  shots.  Also  included  are  data  from  the  1962  North 
Carolina  shots.  The  consistency  of  the  smaller  residuals  at  a  distance  of  400  to  450  km,  obtained 
at  a  station  located  in  the  Piedmont,  is  indicative  of  a  thinner  upper  crust  in  that  region.  The 
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larger  residuals  ai  550  to  650  km  again  reflect  a  thicker  crustal  section  under  the  Appala¬ 
chians. 

Figure  4  is  a  composite  of  the  travel-time  residuals  of  the  ECOOE  5-  and  10- ton  shots 
of  the  northern  and  southern  profiles  and  the  CHASE  HI,  IV,  and  VII  shots.  A  smaller  scale 
is  used  to  show  the  data  obtained  at  distances  between  1000  and  2000  km.  Of  particular  interest 
in  this  figure  is  the  apparent  increase  in  P  velocity  to  8.5  km/sec  at  distances  beyond  1100  km. 

A  similar  velocity  increase  at  comparable  distances  was  noted  by  other  participants  of  the 
EARLY  RISE  program.  This  increase  in  velocity  also  coincides  with  a  substantial  increase  in 
particle- velocity  amplitude  beyond  1100  km.  Across  the  Michigan  Basin  the  first  arrival  ampli¬ 
tudes  were  approximately  10  to  15  dB  larger  than  the  signals  recorded  at  The  University  of 
Michigan's  Botanical  Garden  Station,  located  near  the  southeast  rim  of  the  Basin  (A  s  975  km). 

Figure  5  shows  the  reduced  travel  times  for  the  shots  fired  by  the  USGS  in  Lake  Superior 
during  the  month  of  October  1964.  A  series  of  1-,  2-,  3-,  6-,  and  10- ton  shots  were  fired  at 
several  locations.  A  crustal  velocity  of  6.17  km/sec  and  an  upper-mantle  velocity  of  8.29  km/sec 
were  obtained. 

Figure  6  shows  the  reduced  travel- time  data  for  the  EARLY  RISE  shots  which  were  fired 
in  Lake  Superior  and  were  recorded  along  the  same  line  as  shown  in  the  previous  figures. 

Many  of  the  same  recording  sites  were  reoccupied.  An  upper- mantle  velocity  of  8.0  was  used 
for  these  computations.  A  crustal  velocity  of  6.3  km/sec  verifies  the  results  shown  in  the  figure 
5.  However,  a  clear  7.7-km  lower  crustal  layer  is  indicated.  (This  is  illustrated  in  figure  9.) 
Considerable  scatter  can  be  seen  in  the  Pn  arrivals  beyond  700  km.  This  is  due  primarily  to 
the  high  background  noise  caused  by  increased  cultural  activity  across  southeastern  Michigan, 
Ohio,  and  West  Virginia.  The  evidence  for  an  Appalachian  mountain- root  system  is  not  as  con¬ 
vincing  on  this  profile.  The  mountains  are  1175  to  1360  km  from  the  shot  point.  The  residuals 
on  the  western  half  of  the  mountains  are  smaller  than  those  on  the  eastern  half  and  those  located 
just  to  the  east  of  the  mountains.  However,  the  scatter  of  data  closer  to  the  shot  point  would 
cause  one  to  think  this  anomaly  could  also  be  due  to  scattering.  For  the  shots  recorded  at 
distances  of  1500  to  1650  km,  the  first  arrival  could  not  be  read  accurately.  At  our  most  distant 
station  the  signal  amplitude  was  larger  even  though  the  background  noise  was  high. 

The  vertical- component  recordings  for  four  shots  recorded  at  the  same  station  are  shown 
in  figure  7.  These  four  individual  shots  were  dubbed  onto  one  magnetic  tape,  taking  into  account 
zero  times  of  the  shots  and  slight  differences  in  shot-point  locations.  The  times  were  aligned 
so  that  each  shot  is  within  ±2  msec  at  the  time  of  the  first  arrival.  The  maximum  deviation  over 
the  entire  record  is  ±5  msec.  The  vertical  components  shown  here  were  filtered  on  playback 
with  a  1-  to  4-cps  passband,  which  greatly  improved  the  signal-to-noise  ratio.  The  summed 
vertical,  longitudinal,  and  transverse  components  for  the  four  shots  are  also  shown  on  figure  7. 

At  this  distance  it  can  be  seen  that  most  of  the  arrival  energy  of  the  first  compressional  wave 
is  recorded  on  the  vertical  component.  The  summing  technique  (type  of  array  processing)  gives 
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a  12-dB  improvement  in  signal-to-noise  ratio  of  the  first  arrival.  Such  techniques  allow  a 
more  accurate  determination  of  the  first  arrival  times. 

Figure  8  shows  a  composite  of  a  section  of  the  EARLY  RISE  recordings.  The  details  at 
this  scale  are  difficult  to  see;  however  the  increased  noise  beyond  750  km  is  quite  evident.  One 
feature  particularly  interesting  is  the  small- amplitude  first  arrivals  at  stations  located  at  160 
to  275  km.  This  is  illustrated  more  clearly  in  figure  9.  At  these  distances  the  high-gain  chan¬ 
nels  were  deliberately  used  to  permit  overrecording;  this  provided  better  record  gains  for  the 
low-amplitude  first  arrivals.  (Figure  9  shows  the  first  arrivals  on  the  vertical  component  with 
the  overrecorded  later  phases  omitted.  The  records  are  adjusted  so  that  arrivals  with  a  veloc¬ 
ity  of  7.7  km/sec  are  aligned  vertically  on  this  section.  The  first  three  traces  are  on  the  6.2- 
km/sec  branch  of  the  travel-time  curve.  If  the  record  gains  used  in  the  field  had  been  set  so 
the  later  phases  would  not  be  overrecorded,  this  early  arrival  would  have  been  very  difficult  if 
not  impossible  to  detect.) 

The  travel-time  data  for  this  reversed  profile  are  shown  in  table  I. 

1.3.  ATTENUATION  DATA 

Figure  10  shows  the  first- arrival  vertical  component  of  ground  displacement  versus  dis¬ 
tance  for  the  EARLY  RISE  data.  Although  there  is  considerable  scatter,  an  average  slope 
proportional  to  the  inverse  3/2  power  of  the  distance  fits  these  data  reasonably  well.  On  the 
reversed  profile,  the  ECOOE  shots,  the  1962  North  Carolina  data,  and  the  CHASE  shots  fit  more 
closely  to  a  curve  whose  slope  is  proportional  to  the  inverse  third  power  of  the  distance.  Fur¬ 
thermore,  size  for  size,  the  shots  fired  in  Lake  Superior  were  more  efficient  in  coupling  energy 
into  seismic  waves  than  the  shots  fired  in  the  Atlantic  Ocean.  The  signal  amplitudes  for  the 
Lake  Superior  shots  ranged  from  6  to  18  dB,  or  a  factor  of  2  to  8  larger  for  shots  recorded  at 
comparable  distances.  A  partial  explanation  for  this  phenomenon  may  be  that  the  shots  in  Lake 
Superior  (especially  in  EARLY  RISE)  were  fired  at  an  optimum  depth  for  signal  enhancement  due 
to  reinforcement  by  the  bubble  pulse  (optimum  depth  between  593  and  607  ft).  However,  it  is 
believed  that  Lake  Superior  is  an  anomalous  region  for  the  coupling  of  seismic  energy. 

For  a  clearer  interpretation  of  the  attenuation  data,  isoparticle- velocity  maps  were  pre¬ 
pared  from  the  spectral  data  of  the  first  compressional-wave  arrival  at  each  station.  This 
technique  showed  that  many  of  the  fluctuations  in  spectral  data  persisted  from  shot  to  shot  where 
the  station  remained  fixed  and  from  station  to  station  when  the  shot  point  was  fixed.  These  trends 
are  difficult  to  follow  when  the  individual  spectral  curves  are  examined.  Figure  11  shows  the 
vertical  component  of  the  ground- particle  velocity  for  the  first  three  cycles  of  the  first  com¬ 
pressional-wave  arrival  recorded  on  the  EARLY  RISE  shots.  The  gap  in  data  from  270  to  380 
km  occurs  where  the  profile  crosses  Lake  Michigan.  The  dots  correspond  to  data  points  (the 
values  have  been  omitted  for  sake  of  clarity).  The  particle- velocity  amplitudes  are  contoured 
at  6-dB  (factor  of  2)  intervals.  These  amplitudes  correspond  to  true  ground  motion.  All  instru¬ 
ment  responses  have  been  corrected  out.  The  spectra  of  the  EARLY  RISE  shots  are  strongly 
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TABLE  I.  TRAVEL-TIME  DATA 


ECOOE 

ECOOE 

(Southern  Profile) 

(Northern  Profile) 

EARLY  RISE 

Vhi 

2.0/0.48 

2. 0/1. 7 

4.75/2.9 

Vl/h2 

6.0/28.7 

6.0/25.8 

6.2/22.5 

v2/h3 

7.7/25.0 

v  /h. 
n  4 

8.2/70.0 

8.2 

8.3/21.6 

V 

8.5 

8.4 

Total  crustal 
thickness 

29.2  km 

27.5  km 

50.4  km 

Depth  to 
high-velocity 
upper -mantle 
layer 

99.2 

72  km 

influenced  by  the  bubble- pulse  frequency  and  the  reverberation  of  the  sound  waves  between  the 
lake  bottom  and  the  surface.  The  bubble-pulse  frequency  is  approximately  2  cps  and  would  be 
in  phase  at  that  water  depth.  The  organ-pipe  frequencies  for  constructive  and  destructive 
interference  are: 

Constructive  Destructive 

fQ  =  2  cps  2fg  =  4  cps 

3fg  =  6  cps  4fjj  =  8  cps 

5fg  =  10  cps  6fjj  =  12  cps 

7fg  =  14  cps 

A  rather  remarkable  correlation  with  these  frequencies  can  be  seen  on  the  particle- velocity 
highs  and  lows  shown  on  figure  11.  The  major  peaks  along  the  entire  profile  occur  at  approxi¬ 
mately  2  and  5-6  cps.  The  former  coincides  with  the  bubble-pulse  frequency  and  the  funda¬ 
mental  organ-pipe  mode,  and  the  latter  coincides  with  the  third  harmonic.  A  high  can  also  be 
correlated  with  the  fifth  harmonic.  The  particle- velocity  lows  correlate  quite  well  with  the  second 
and  fourth  harmonic,  where  destructive  interference  should  occur.  A  particle-velocity  low  at  all 
frequencies  can  be  seen  near  the  crossover  point  (160  km)  of  the  6.2-  and  the  7.7-km/sec  branch¬ 
es  of  the  travel-time  curves.  Unfortunately,  the  Pn  cross  overpoint  falls  near  that  section  of 
Lake  Michigan  where  we  have  no  data.  One  other  predominant  feature  should  be  pointed  out.  A 
pronounced  particle- velocity  high  occurs  at  approximately  885  km.  This  same  high  was  ob¬ 
served  at  a  group  of  four  stations  located  near  Anna,  Ohio,  at  the  same  distance.  The  latter 
area  is  150  km  southwest  of  the  main  profile  and  Is  an  area  where  small  earthquakes  have  been 
recorded  in  the  past.  The  cause  of  this  large  amplitude  is  not  known.  It  may  possibly  corre¬ 
late  with  a  cusp  in  the  travel-time  curves.  Its  buildup  is  observable  at  seven  stations. 

Figure  12  shows  the  isoparticle-velocity  map  for  the  1963  Lake  Superior  south  line.  This 
profile  was  obtained  by  two  crews  which  changed  sites  each  night.  They  recorded  nine  different 
shots  along  a  line  perpendicular  to  the  center  of  the  line  of  shot  points.  The  particle- velocity 
highs  correlated  well  with  the  bubble-pulse  and  organ-pipe  modes  that  would  cause  reinforce¬ 
ment.  The  predominant  signal  along  the  entire  profile  is  4  cps,  which  correlates  with  the  first 
odd  harmonic  of  the  organ-pipe  mode.  The  sharp  decrease  in  particle- velocity  amplitude  at 
160  to  170  km  correlates  with  the  data  on  the  EARLY  RISE  shot  made  at  the  crossover  point  on 
the  travel-time  curve  from  the  6.2  to  7.7  branch. 


In  order  to  investigate  the  attenuation  of  a  seismic  wave,  using  the  type  of  data  shown  in 
this  figure,  an  equation  of  the  type 


A  -A  X-e'^ 
Aij  "  ojXi 


(1) 


was  employed.  It  accounts  for  geometrical  spreading  (n)  and  for  the  absorption  (a)  of  the  seis¬ 
mic  waves.  A^  represents  the  particle-velocity  amplitude  for  the  i-th  shot  and  the  j-th  frequency; 
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Aq.  represents  a  constant  assumed  to  be  independent  of  distance  but  dependent  on  the  j-th  fre¬ 
quency;  X.  represents  the  distance  from  the  seismometer  to  the  i-th  shot;  and  f^  is  the  j-th 
frequency. 

To  derive  equation  1,  we  begin  with 


-a  f;X. 


A..  =  A  .X.  “  e 

l]  l]  l 


and  solve  for  A. . : 

i] 


a  f.X. 

v-..,.x"e  *  > 

ij  i]  i 


A0.  =A 


fi  =  1,  2, 

li  =  1,  2, 


If  we  now  let 


A0. 

i] 


M 

N 


where  M  is  the  number  of  shots,  we  have  arrived  at  equation  1.  Hence,  for  each  A^.,  an  Ay  was 
calculated  so  that  a  matrix  was  obtained  whose  rows  (1,  2,  .  .  .  ,  j,  .  .  .  ,  N)  represent  the  N  dif¬ 
ferent  frequencies,  and  whose  columns  (1,  2,  .  .  .,  i,  .  .  .,  M)  represent  the  M  different  shot 
distances  (in  increasing  order).  An  average  Aq,  was  found  for  each  frequency  across  the  dis¬ 
tances,  by  summing  the  A°  for  a  constant  j  and  dividing  the  sums  by  M.  There  are  then  N  Aq/s, 
which  are  in  a  respect  independent  of  distance.  This  matrix  represents  an  average  spectrum 

<Ao}»- 

Using  equation  1  and  the  computed  Ao.,  new  A^'s  (designated  A^)  are  computed,  which  in 
general  will  be  different  from  the  original  A.  .'s. 

The  Ajj's  and  Aj.'s  are  converted  from  particle- velocity  amplitudes  to  the  equivalent  ampli¬ 
tudes  in  decibels  by  using  the  following  relationship: 


Aij  (in./sec)  =  20  log^-gL^- 

A  matrix  of  residuals  (Sj.)  are  formed  where 

S..  =  A..  -  a:. 

i]  ij  i] 


The  Sy's  are  plotted  on  a  graph  whose  abscissa  is  the  distance  and  whose  ordinate  is  the 


frequency. 


Figure  13  shows  the  summation  of  the  residuals  versus  the  absorption  coefficient  for  values 
of  n,  the  geometrical  spreading  factor,  ranging  from  0  to  3.  The  best  fit  for  n  and  a  would  cor- 

_4 

respond  to  the  residual  minimum.  For  this  profile  a  value  of  n  =  2.5  and  a  =  7  x  10  yields  the 
smallest  residual.  However,  a  value  of  n  =  1.5  would  correspond  to  geometrical  spreading  due 

-3 

to  head-wave  propagation  (perhaps  a  value  of  n  =  2  for  these  distances).  For  n  =  1.5,  a  =  1  x  10 
gives  the  best  fit.  9 


A 


In  the  foregoing  discussion  we  have  assumed  that  the  absorption  would  be  proportional  to 
the  first  power  of  the  frequency.  In  order  to  examine  this  assumption,  residuals  were  computed 

for  I**,  where  P  ranged  from  0  to  2.  The  results  are  shown  in  figure  14.  It  can  be  seen  that 

1/2  -3 

f  gives  the  smallest  residual  at  a  =  5  x  10  for  n  =  1.0. 

Figure  15  shows  the  reduced  isoparticle- velocity  residuals  after  the  geometrical  spreading 

.  3 

factor  and  absorption  factor  have  been  removed.  In  this  case  values  of  n  =  1.0  and  a  =  1.5  x  10 
were  used.  This  does  not  give  the  best  fit  from  the  data  shown  in  figure  14,  but  it  does  show  the 
technique  that  can  be  used  to  help  study  the  anomalous  amplitudes  encountered  along  the  profile. 
The  bubble-pulse  frequency  is  indicated  by  the  X's  on  this  graph. 

Figure  16  shows  the  particle- velocity  amplitudes  for  the  150- km  site  on  the  ECOOE  southern 
profile.  Most  of  those  data  fall  on  the  6.0  branch  of  the  travel- time  curve.  The  relief  or  change 
in  signal  here  is  less  than  that  of  the  EARLY  RISE  or  1963  Lake  Superior  data.  One  reason  is 
that  the  bubble- pulse  frequency  is  below  1.5  cps  for  all  but  the  two  furthest  shots  and  the  funda¬ 
mental  organ-pipe  frequency  is  above  10  cps  for  all  but  the  last  three  shots.  The  particle-veloc¬ 
ity  highs  for  the  last  three  shots  correlate  with  the  bubble- pulse  frequency  and  the  fundamental¬ 
mode  organ- pipe  frequency. 

Figure  17  shows  the  sum  of  the  residuals  versus  a  for  n  ranging  from  0.5  to  3  for  the  ECOOE 

_4 

southern- profile  150-km  site.  The  best  fit  is  for  n  =  0.5  and  a  =  6  x  10  .  For  n  =  1.5,  «  =  1  x 

.4 

10  would  give  the  best  fit. 

From  data  of  this  type  it  is  anticipated  that  we  will  be  able  to  determine  the  mathematical 
expression  that  will  be  best  suited  to  predict  signal  amplitudes  versus  distance  and  to  make 
accurate  determinations  of  crustal  and  mantle  Q's. 
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FIGURE  1.  MAP  OF  SHOT  POINTS  AND  RECORDING  STATIONS 


»E  SERIES  AND  CHASE  ID,  IV,  AND  VH 


FIGURE  5.  TRAVEL-TIME  RESIDUALS  FOR 
OCTOBER  1964  LAKE  SUPERIOR  EXPERI¬ 
MENT 
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FIGURE  6.  TRAVEL-TIME  RESIDUALS  FOR  THE  1966  EARLY 
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FIGURE  10.  GROUND  DISPLACEMENT  VS.  DISTANCE  FOR  PROJECT  EARLY  RISE 
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FIGURE  12. 


ISOPARTICLE -VELOCITY  MAP  OF  THE  FIRST  P-WAVE  ARRIVAL  FOR  1963 
LAKE  SUPERIOR  SOUTH  LINE 
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FIGURE  16.  ISOPARTICLE-VELOCITY  MAP  FOR  FIRST  P-WAVE  ARRIVALS  ON  ECOOE 
150-KM  STATION  ON  SOUTHERN  PROFILE 
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THEORETICAL  STUDIES  IN  ELASTIC-WAVE  SCATTERING 
I.  K.  Mclvor 


To  evaluate  the  effect  of  nonhomogeneities  on  the  propagation  of  elastic  waves,  a  contin¬ 
uing  analytical  investigation  of  elastic-wave  scattering  is  being  carried  out.  The  initial  investi¬ 
gation  considered  the  scattering  of  plane  waves  in  an  elastic  media  by  elliptical  obstacles  or 
cavities.  A  general  solution  was  obtained,  and  illustrative  numerical  calculations  were  made 
for  the  degenerate  case  of  a  slit.  It  was  found  that  the  scattered  field  exhibited  strong  direc¬ 
tional  selectivity. 

The  major  area  of  investigation  at  present  is  the  scattering  effect  of  small  imperfections 
at  the  interface  of  elastic  media.  A  perturbation  method  has  been  developed  which  is  applicable 
to  a  wide  range  of  problems.  A  specific  solution  has  been  obtained  for  the  scattering  of  body 
waves  by  surface  irregularities  on  a  half  space.  The  solution  is  expressed  as  a  surface  convo¬ 
lution  of  a  known  kernel  with  the  imperfection.  Thus  numerical  results  for  arbitrary  imperfec¬ 
tion  shapes  are  readily  obtained.  In  the  present  problem  the  surface  displacements  and  total 
scattered  energy  have  been  computed  for  a  number  of  different  scatterers. 

It  is  anticipated  that  the  method  will  be  applied  to  other  problems  of  seismic  interest, 
including  the  scattering  of  surface  waves  and  the  effect  of  imperfections  on  the  reflection  and 
transmission  of  body  waves  at  a  material  discontinuity.  Finally,  an  effort  needs  to  be  made  to 
compare  the  theoretical  results  with  actual  seismic  data  to  see  if  major  structural  imperfections 
can  account  for  observed  anomalies. 
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LAKE- BOTTOM  SEISMIC-NOISE  STUDIES 
D.  E.  Willis 

The  series  of  underwater  shot  programs  conducted  in  Lake  Superior  during  1963,  1964, 
and  1966  rangtd  in  size  from  1  to  10  tons  of  high  explosives  and  have  been  detected  at  distances 
in  excess  of  2000  km.  Underwater  shots  provide  a  very  efficient  coupling  mechanism  for  seis¬ 
mic  waves.  The  increased  signal  level  for  a  properly  tamped  water  shot  can  generate  ground 
displacements  that  are  a  factor  of  20  or  more  larger  than  equivalent  shots  fired  in  rock.  The 
seismic  signals  generated  by  these  shots  in  Lake  Superior  were  found  to  have  larger  amplitudes 
(up  to  a  factor  of  8)  at  comparable  distances  than  equivalent  shots  fired  in  the  Atlantic  or  the 
Pacific  Oceans. 

Since  the  Lake  Superior  area  was  found  to  have  an  efficient  coupling  of  seismic  energy  it 
was  deemed  highly  desirable  to  ascertain  whether  or  not  there  would  be  any  reciprocity  in  the 
recording  of  teleseismic  events;  that  is,  whether  seismic  sources  generated  at  teleseismic 
distances  in  other  areas  would  cause  larger  ground  displacements  in  the  vicinity  of  the  shot 
points  in  Lake  Superior  where  particularly  good  coupling  was  observed.  The  EARLY  RISE 
shot  point,  located  near  shot  point  51  (see  fig.  18)  of  the  1963  Lake  Superior  Experiment,  was 
chosen  as  the  site  to  conduct  these  tests. 

Figure  18  is  a  map  of  the  portion  of  Lake  Superior  that  is  of  interest  to  this  study.  The 
shot  points  of  the  1963  experiment  are  shown  along  with  the  water  depths  and  contours  of  the 
water  depth  [2],  The  contours  are  not  in  complete  agreement  with  the  indicated  water  depths. 
The  Obenhoff  land  station  located  on  the  Keweenaw  Peninsula  is  also  indicated.  The  dashed 
line  shown  in  this  figure  is  the  location  of  a  sub-bottom  seismic  reflection  profile  made  by  Dr. 

R.  J.  Wold  of  the  University  of  Wisconsin.  This  profile  is  shown  in  figure  19. 

A  pronounced  structural  high  having  approximately  350  ft  of  relief  was  disclosed  by  this 
profile.  This  feature  is  approximately  9  km  in  diameter  and  has  a  more  or  less  flat  plateau  on 
its  southeast  flank.  This  plateau  was  chosen  as  the  site  to  install  a  lake-bottom  seismograph 
since  the  lake- bottom  sediments  appeared  to  be  thin.  The  plateau  has  a  water  depth  of  approxi¬ 
mately  530  to  550  ft. 

The  lake- bottom  seismograph  was  obtained  on  loan  from  Professor  Hugh  Bradner  of  the 
Institute  of  Geophysics  and  Planetary  Physics,  University  of  California  at  San  Diego,  LaJolla. 

A  detailed  description  of  this  unit  may  be  found  in  reference  3.  The  package  consists  of  matched 
three- component  equiangular  seismometers,  amplifiers,  power  supplies,  a  tape  recorder,  and 
an  electronic  clock.  The  following  modifications  of  these  units  were  made  to  adapt  them  to  the 
purposes  of  this  investigation.  The  high-frequency  response  of  the  amplifiers  and  tape  recorder 
was  extended  to  17  cps;  the  seismometers  were  tuned  to  a  natural  period  of  1.1  sec  with  0.7 

28 


critical  damping;  the  power  supply  was  modified  to  provide  13  days  of  recording  time;  the 
interval  clock  was  modified  to  provide  1-min  and  1-hr  time  pulses  with  a  10-cps  carrier  last¬ 
ing  1  sec  and  1  min,  respectively. 

The  response  of  the  seismometers  used  in  the  lake- bottom  unit  and  the  Obenhoff  land  station 
is  shown  in  figure  20.  The  land-station  seismometers  (matched  three- component)  have  approxi¬ 
mately  14  dB  more  output  at  1  to  2  cps  and  22  dB  at  10  cps  than  the  lake- bottom  unit. 

The  lake- bottom  seismograph  was  tethered  in  the  manner  shown  in  figure  21.  Nylon  line 
was  used  to  secure  the  seismograph,  anchor,  and  buoy.  Ample  slack  was  used  to  reduce  the 
buoy  motion  being  transmitted  through  the  nylon  line  to  the  seismograph. 

The  package  was  emplaced  with  the  assistance  of  the  USCG  cutter  WOODRUSH  and  was 
operated  for  three  periods  during  June- August  1967.  The  times  of  operation  are  shown  in 
table  II. 

Very  high  background  noise  was  encountered  on  all  three  drops.  Long  periods  of  high 
background  noise  were  found  that  can  be  correlated  with  meteorological  conditions  and  with 
wave  action.  This  is  illustrated  in  table  III  and  figures  22  and  23.  Table  III  shows  the  average 
seismometer  output  voltage  from  the  three- component  lake- bottom  seismograph  and  the  three- 
component  land-station  seismograph  at  4  hr  intervals  for  the  time  period  0800,  June  7,  to  0800, 
June  10,  1967,  Also  shown  is  the  wind  velocity  and  direction  as  well  as  the  lake-wave  height 
and  direction.  No  correction  has  been  made  for  the  difference  in  seismometer  response  for  the 
land  and  lake-bottom  seismographs,  so  this  difference  should  be  kept  in  mind  when  responses 
at  the  two  stations  are  compared. 

The  periods  of  peak  wind  velocity  and  lake-wave  height  which  occur  between  0800  and  1600 
on  June  8  correlate  quite  closely  with  the  maximum  signal  levels  at  both  the  land  station  and  the 
lake-bottom  unit.  The  three  lake-bottom  seismometers  have  nearly  identical  outputs  at  any 
given  time.  The  land- station  seismograph  displays  the  lowest  levels  on  the  vertical  component 
and  the  highest  on  the  east-west  component.  The  prevailing  winds  could  account  for  the  latter. 

Figure  22  shows  the  variation  in  average  background  noise  at  15-min  intervals  for  the  time 
period  1000,  June  7,  to  2330,  June  8.  The  land- station  vertical- component  seismometer  output 
is  shown.  Since  the  output  of  the  land-station  seismometers  is  14  to  22  dB  greater  than  that  of 
the  lake- bottom  seismometers,  it  can  be  seen  that  the  difference  in  noise  levels  between  these 
two  sites  is  quite  large. 

This  slow  buildup  and  decay  of  background  noise  was  found  to  be  rather  common  during  the 
three  drop  periods.  The  maximum  levels  persisted  for  intervals  ranging  from  several  hours 
to  several  days. 

Spectral  analysis  of  the  average  background  noise  for  four  1-hr  intervals  is  shown  in  figure 
23.  These  time  intervals  include  the  quiet  period  prior  to  the  buildup,  the  maximum  noise  levels, 
and  the  intermediate  levels  on  buildup  and  decay.  Several  significant  features  can  be  seen  in 
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TABLE  III.  DATA  I  OR  THE  THREE  DROPS 
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this  figure.  The  minimum  and  maximum  velocities  at  1  cps  are  indicated  for  both  stations. 

The  peak  ground  displacement  ranges  from  33  to  1500  mji  on  the  lake  bottom  and  from  0.4  to 
10  mp  on  the  land  station.  These  stations  are  located  approximately  60  km  apart.  It  can  be 
seen  that  the  lake- bottom  noise  is  nearly  two  orders  of  magnitude  higher  than  that  at  the  land 
station.  The  lake- bottom  noise  has  three  prominent  peaks  occurring  at  0.5,  1.0,  and  5.0  cps, 
with  a  pronounced  minimum  at  2  to  3  cps.  The  land  station  has  a  pronounced  peak  at  2.5  cps 
with  a  minor  peak  at  5  cps  and  a  pronounced  minimum  at  1.6  to  2  cps.  The  peak  on  the  land 
station  is  centered  near  the  minimum  on  the  lake- bottom  station. 

Figure  24  shows  a  similar  buildup  and  decay  in  background  noise  on  the  third  drop.  The 
shapes  of  the  spectral  curves  are  similar  to  those  of  figure  23,  with  the  higher  frequency  peak 
occurring  at  4  cps.  A  second  lake- bottom  seismograph  was  also  used  on  the  third  drop.  This 
unit  was  a  free-fall  untethered  station  that  was  used  to  compare  background  no;se  with  the 
tethered  unit.  No  significant  difference  in  sif  ial  level  was  observed. 

Shown  in  figures  25  and  26  are  seismograms  obtained  with  both  the  Lake  Superior  units 
during  the  Buffalo,  New  York,  earthquake  (magnitude  3.9)  of  June  13,  1967.  The  bottom  trace  of 
figure  25  is  the  summed  output  of  the  three  lake-bottom  seismometers  and  would  correspond 
to  a  true  vertical  component.  Two  significant  facts  can  be  observed  on  lese  seismograms. 
Taking  into  account  the  difference  in  response  of  the  seismometer  and  the  difference  in  record 
gains,  the  ground  amplitudes  on  the  lake- bottom  unit  are  larger  than  those  on  the  land- station 
seismograph.  The  time  duration  of  the  event  on  the  lake- bottom  unit  is  a  factor  of  2  or  more 
longer  than  that  of  the  land  station.  Spectral  analyses  of  this  event  show  a  pronounced  peak  at 
1  cps.  The  important  factor  in  comparing  these  two  recordings  is  the  signal- to- noise  ratio  as 
a  function  of  frequency.  This  is  shown  in  figure  27.  It  can  be  seen  from  this  figure  that,  although 
the  signal  level  on  the  land  station  was  lower  than  that  on  the  lake- bottom  station,  the  maximum 
signal-to-noise  ratio  was  obtained  by  the  land  station. 

For  comparison  purposes  the  signal-to-noise  ratios  for  this  same  event  recorded  at  The 
University  of  Michigan's  Botanical  Garden  intermediate- depth  well  station  and  at  the  surface  in 
Ann  Arbor,  Michigan,  are  also  shown  in  figure  27.  The  superiority  of  the  well  recording  is  quite 
evident  in  spite  of  the  fact  that  the  well  is  not  of  the  proper  depth  (the  surface  reflection  inter¬ 
ference  for  the  seismometers  at  a  depth  of  1206  feet  causes  a  minimum  at  2-2.5  cps). 

The  Denver  earthquake  (magnitude  5.3)  of  August  9,  1967,  was  also  recorded  by  the  lake- 
bottom  and  the  land  seismograph  units.  The  spectral  curves  and  signal-to-noise  ratios  were 
similar  to  those  of  the  Buffalo  earthquake.  Pronounced  peaks  at  1  cps  were  observed  on  the 
lake-bottom  unit. 

Numerous  quarry  shots  were  also  recorded  by  both  units.  A  striking  feature  of  these  local 
events  was  a  very  sharp  spectral  peak  at  4  cps  on  the  lake- bottom  unit.  The  seismograms  dis¬ 
play  a  pronounced  ringing  effect  that  were  very  dissimilar  to  the  land-station  recordings. 
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The  lake -bottom  units  and  the  land-station  equipment  have  been  used  to  record  teleseismic 
and  local  events  at  the  Geophysics  Laboratory  (near  Ypsilanti,  Michigan).  Absolute  signal  levels 
and  the  spectral  data  from  both  units  are  in  very  good  agreement.  Hence,  if  the  lake- bottom 
seismograph  "plant"  is  not  introducing  a  resonance  effect,  the  pronounced  ringing  and  sharp 
spectral  peaks  for  local  events  probably  can  be  correlated  with  some  sort  of  water  reverberation 
effect.  Teleseismic  events  with  peak  signals  below  4  cps  do  not  produce  this  sharp  ringing  effect. 
However,  the  distant  events  recorded  on  the  lake- bottom  unit  persist  for  longer  time  intervals 
than  on  the  land  station. 

It  can  be  concluded  that  the  operation  of  lake  bottom  seismographs  in  Lake  Superior  would 
not  be  desirable  on  a  permanent  basis  because  of  the  poor  signal-to- noise  ratios  compared  with 
those  of  land  stations,  because  of  long  periods  of  very  high  background  noise,  and  because  of  the 
ringing  character  of  the  seismograms. 
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FIGURE  18.  LOCATION  MAP  OF  LAKE-BOTTOM  AND  LAND-STATION  SEISMO¬ 
GRAPHS,  LAKE  SUPERIOR 
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FIGURE  20.  LAKE  BOTTOM  AND  LAND  SEISMOMETER  RESPONSE  CURVES,  LAKE  SUPERIOR 


FIGURE  22.  COMPARISON  OF  SEISMIC  BACKGROUND  NOISE  AT  15-min  INTERVALS  FOR 
JUNE  7  AND  8,  1967,  LAKE  SUPERIOR  LAND  STATION 
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FIGURE  23.  COMPARISON  OF  SEISMIC  BACKGROUND  NOISE  FOR  FOUR  INTERVALS, 

BOTH  LAKE  SUPERIOR  STATIONS 
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FIGURE  27.  SIGNAL-TO-NOISE  RATIO  OF  BUFFALO,  NEW  YORK,  EARTHQUAKE  OF 

JUNE  13,  1967 
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SEISMIC  RADIATION  FROM  SHOAL 
C.  G.  Bufe  and  D.  E.  Willis 

Travel  times,  amplitude,  and  particle- velocity  spectra  for  the  SHOAL  event  of  26  October 
1963  have  been  examined  for  distances  to  1200  km  and  for  widely  distributed  azimuths.  A 
composite  curve  of  raw  travel-time  data  reaffirms  the  low  Pn  velocities  previously  observed 
in  the  Basin  and  Range  province  and  yields  a  crustal  thiciTiess  of  about  26  km.  The  travel- 
time  data  was  obtained  from  LRSM,  VELA  array,  USCG,  Sandia  Corporation,  University  of 
California,  University  of  Michigan,  and  various  volunteer  team  sources.  These  data,  when 
corrected  for  station  elevations,  yield  a  least- square  Pn  velocity  of  7.84  km/sec.  Residuals  to 
this  Pn  curve  have  been  contoured  for  distances  from  200  to  1200  km  and  are  shown  in  figure 
28.  With  the  exception  of  the  central  California  and  southern  Nevada  data,  these  residuals  bear 
a  negative  correlation  to  the  mean  station  residuals  obtained  by  Carder,  Gordon,  and  Jordan 
[4]  for  explosive  sources  at  teleseismic  distances,  and  thus  seem  to  reflect  lateral  variations 
in  mantle  velocity  just  below  the  crust.  Paths  through  Herrin  and  Taggart's  [5]  high  and  low 
Pn  apparent  velocity  zones  explain  the  observed  residual  pattern  north  and  east  of  SHOAL. 

Amplitude  data  from  LRSM,  VELA  array,  Sandia,  and  University  of  Michigan  stations  were 
used  in  a  study  of  the  attenuation  of  Pn  at  distances  of  200  to  1200  km.  A/T  was  found  to  decrease 
as  the  inverse  cube  of  the  distance  (A  =  amplitude;  T  =  period).  Residuals  from  the  mean  in¬ 
verse  cube  were  measured  in  decibels  and  are  shown  in  figure  29.  Negative  residuals  are  within 
the  shaded  areas.  There  is  a  zone  of  large  positive  Pn  amplitude  residuals  in  northeastern 
Nevada,  but  generally  the  first  arrivals  to  the  east  of  SHOAL  are  lower  in  amplitude  than  those 
at  the  same  distances  to  the  south  and  north.  Amplitude  data  from  SHOAL  and  from  the  Fallon 
earthquake  of  July  20,  1962,  contoured  by  Jordan,  Black,  and  Bates  [6]  also  indicate  a  pattern 
of  anomalously  high  Pn  amplitudes  in  northeastern  Nevada  and  northern  Utah. 

Two  sets  of  three  stations  to  the  north,  northeast,  and  east  of  SHOAL  at  distances  of  320  and 
1100  km  were  chosen  for  analysis  of  particle-velocity  spectra.  Four  of  these  were  portable 
stations  operated  by  The  University  of  Michigan  and  two  were  LRSM  vans  operated  by  the  Geo¬ 
technical  Corporation. 

The  spectra  of  Pn  at  The  University  of  Michigan's  320- km  stations  are  shown  in  figure  30. 
Dinner  Station,  Nevada,  which  is  northeast  of  SHOAL,  displays  much  more  energy  below  2.5  Hz 
than  the  stations  to  the  north  or  east.  Above  3  cps,  Chicken  Springs,  Oregon  (north  of  SHOAL), 
and  Dinner  Station  have  particle  velocities  well  above  those  at  Berry  Creek  to  the  east.  The 
maximum  P-particle  velocities  (fig.  31)  show  a  stronger  similarity  between  stations,  each  peak¬ 
ing  at  2.5  cps,  although  Chicken  Springs  sees  more  signal  above  3.0  cps  than  the  other  320-km 
stations. 
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The  Pn  spectra  at  distances  near  1100  km  are  shown  in  figure  32.  At  all  frequencies 
shown,  the  Tres  Piedres  site  to  the  southeast  shows  lower  amplitudes  than  the  other  two  sta¬ 
tions.  The  spectra  of  maximum  P  as  shown  in  figure  33  are  similar  at  the  three  stations,  al¬ 
though  Shepherd  displays  a  lower  signal  level  between  1  and  2  cps  than  the  other  stations.  Van 
Nostrand  and  Helterbran  [7]  have  used  data  from  other  stations  to  compute  the  power  spectra 
for  the  SHOAL,  HARDHAT,  and  "Fallon  earthquake"  events. 
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MODE  FILTERING 
R.  M.  Turpening 

With  the  optical  data  processing  system  as  a  guide,  work  has  progressed  toward  the  under¬ 
standing  of  present  linear  mode  filters,  development  of  new  ones,  and  digital  application  for  all 
of  them.  Much  of  the  work  was  performed  while  the  author  participated  in  a  postdoctoral  pro¬ 
gram  at  Massachusetts  Institute  of  Technology,  and  that  opportunity  is  gratefully  acknowledged. 

The  geometry  of  the  system  is  the  same  as  that  given  in  reference  8.  Detailed  digital  anal¬ 
ysis  of  the  square  of  the  two-dimensional  amplitude  "spectrum"  (diffraction  patterns)  of  synthetic 
Rayleigh  waves  (noise)  and  P  waves  shows  that  fine  structure  exists  that  was  not  seen  in  optical 
processing.  Figures  34,  35,  and  38  indicate  that  there  exist  useful  differences  between  signal 
and  noise  "spectra"  at  least  down  to  a  signal- to- noise  ratio  of  1/2  on  the  vertical  seismometer. 

The  contours  are  logarithmic,  and  the  dashed  contours  are  areas  where  the  noise  spectrum 
exceeds  the  signal  spectrum. 

The  filter  which  seems  to  be  indicated  by  the  spectral  differences  results  ir.  the  following 
expression  for  a  filter  output: 

P(a,  t)  =  A  (a,  /3)Z(t)  +  B(a,  0)Z(t)  +  C  (a,  /3)DN(t)  +  D(a,  0)DN(t)  +  E(a,  p)DTN(t)  +  F(a,  0)DTN(t) 

where  Z(t)  is  vertical  time  trace 

DN(t)  is  time  trace  in  direction  of  noise 
DTN(t)  is  time  trace  in  direction  of  transverse  to  the  noise 

A,  B,  C,  D,  E,  and  F  are  weighting  functions  derived  from  the  location  of  the  filter  bands 

in  the  w-(3  plane  and  the  geometry. 

~  is  an  indication  of  a  constant  phase  shift  of  jt/2 

a  is  index  variable  for  computed  projections 

0  is  transformed  variable  related  to  a 

u  is  frequency  in  radians  per  second 

Since  I  (t)  does  not  represent  a  phase  shift  of  n/2  for  all  frequencies,  phase  distortion  will 
result.  Therefore,  it  is  envisioned  that  the  banded  filter  structure,  shown  in  figures  34,  35,  and 
36,  be  applied  repeatedly  along  some  path  through  the  w-(3  plane.  Figure  37  represents  the  two- 
dimensional  spectrum  (diffraction  pattern)  of  a  typical  seismic  event  and  thus  describes  a  possi¬ 
ble  path  for  the  repeated  use  of  the  banded  filter  structure. 


52 


FIGURE  34.  THEORETICAL  DIFFERENCE  OF  SIGNAL  AND  NOISE  TWO- 
DIMENSIONAL  AMPLITUDE  "SPECTRA"  (DIFFRACTION  PATTERNS)  FOR 
NOISE-ELLIPSE  HALF  AXES  EACH  OF  SIZE  1.  Signal  amplitude  is  2  and 
angle  of  incidence  is  45°.  Suggested  filter  is  shown. 
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FIGURE  25.  THEORETICAL  DIFFERENCE  OF  SIGNAL  AND  NOISE  TWO- 
DIMENSIONAL  AMPLITUDE  "SPECTRA”  (DIFFRACTION  PATTERNS)  FOR 
NOISE -ELLIPSE  HALF  AXES  OF  SIZES  3/2  AND  1.  Signal  parameters 
same  as  figure  34.  Suggested  filter  is  shown. 
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S^;»T»E.?RETICAL  difference  of  signal  and  noise  two- 
dimensional  AMPLITUDE  "SPECTRA"  (DIFFRACTION  PATTERNS)  FOR 
NOISE-ELLIPSE  HALF  AXES  OF  SIZES  2  AND  1.  sSfpLa  Jelers  iamf 
as  figure  34.  Suggested  filter  is  shown. 
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FIGURE  37.  DIFFRACTION  PATTERN  OBTAINED  FROM  A  TYPICAL  SEISMIC  EVENT 
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LASA  DATA  CONVERTOR 
P.  L.  Jackson 

A  digital  convertor  for  LASA  data  has  been  completed  and  is  now  operational.  The  convertor 
is  contained  in  one  7- ft  rack  except  for  the  digital  tape  playback  unit.  Fifty  data  channels  can 
be  selected  from  the  LASA  digital  tape.  This  limited  number  of  data  channels  was  selected  for 
reasons  of  economy,  and  can  be  augmented  to  as  many  as  all  the  data  channels,  with  additional 
expense. 

The  convertor  is  a  flexible  research  tool  in  which  LASA  data  output  can  !  scanned  or 
preselected,  viewed  either  one  dimensionally  or  two  dimensionally  in  wiggly  line  or  variable 
density,  and  photographed  as  frame-by-frame  motion  pictures,  continuous  moving  strips,  or 
still  film.  The  output  is  currently  designed  for  a  CRT.  Hard  copy,  e.g.,  50  simultaneous  wiggly 
line  channels,  is  immediately  obtainable  with  Polaroid  photography  of  the  CRT.  As  the  convertor 
can  be  manually  controlled  to  select  and  replay  any  portion  of  the  LASA  tape  and  visually  present 
selected  outputs  for  the  research  seismologist,  it  enables  direct,  interactive,  contact  with 
LASA  data.  He  has  immediately  available  both  dynamic  and  static  presentations  in  several 
different  forms.  Thus  the  convertor  represents  a  powerful  selective  and  interpretive  tool.  When 
further  machine  computation  is  desired,  the  selected  data  segments  can  be  used  with  the  Uni¬ 
versity's  IBM  360/67  computer. 

As  described  in  references  1  and  8,  an  investigation  of  optical  techniques  for  initial  on-line 
detection  of  significant  events  was  planned  upon  completion  of  the  convertor.  However,  the 
decision  to  use  an  entire  IBM  model  40  in  addition  to  a  model  2250  for  detection  alone  will  re¬ 
quire  reassessment  of  this  plan.  Although  optical  techniques,  if  proved  applicable,  would  cost 
a  fraction  of  that  of  model  40,  the  versatility  and  accuracy  would  not  appear  to  be  comparable. 

Figures  38  and  39  show  the  convertor  and  a  closeup  of  selection  panel.  Figure  40  shows 
one  form  of  output. 
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SEISMIC- RAY  TRACING  AND  TRAVEL  TIMES 
BY  DIGITAL  MODELING 
P.  L.  Jackson 

7.1.  INTRODUCTION 

The  purpose  of  this  work  is  to  produce  seismic-ray  traces  through  any  media,  no  matter 
how  inhomogeneous,  with  output  in  the  form  of  travel  times  and  amplitudes.  The  program  under 
developme.it  is  for  both  two-  and  three-dimensional  representations  of  media.  Need  for  such 
programs  has  recently  become  more  severe  because  more  precise  seismology  has  indicated 
lateral  inhomogeneities  in  almost  every  section  and  depth  of  the  earth.  The  program  under 
development  will  enable  us  to  determine  travel  times  for  the  entire  earth  or  for  any  portion  or 
depth  of  the  earth.  For  example,  a  simple  layered  structure  can  be  easily  constructed  and 
later  perturbed  in  selected  areas  to  more  accurately  reproduce  observed  travel  times. 

The  investigation  is  to  provide  the  seismologist  with  a  tool  to  experimentally  determine 
velocity  and  density  structures.  It  was  thought  that  the  digital  computer,  with  its  high  speed, 
versatility,  ease  of  varying  both  the  analysis  and  the  data  inputs,  and  particularly  the  capability 
of  realistic  display  of  its  output  with  a  CRT,  offered  a  unique  potential  for  this  purpose.  The 
aim  was  to  develop  a  means  of  ray  tracing  through  which  travel  times  could  be  computed  and 
displayed,  in  which  the  rays  would  be  displayed  along  with  the  velocity  structure.  With  such 
a  display  and  the  capacity  to  immediately  retrace  the  rays  after  alteration  of  the  indicated 
structure,  an  interactive  tool  for  the  interpretative  seismologist  would  be  achieved. 

Until  recently,  ray  tracing  has  been  confined  to  a  spherically  symmetric  earth,  horizontal 
layers,  or,  in  optics,  through  layers  which  can  be  represented  as  conic  sections.  Spherically 
symmetric  ray  tracing  has  been  described  by  Julian  and  Anderson  [9]  and  shown  by  Lewis  and 
Meyer  [10].  Iyer  and  Punton  [11]  have  described  a  method  to  compute  rays  through  a  two- 
dimensional  inhomogeneous  velocity  field  first  by  the  construction  of  wave  fronts.  Yacoub  [12] 
has  traced  rays  through  inhomogeneous  media  and  has  computed  amplitude  by  means  of 
Zoeppritz's  equations.  However,  information  on  Yacoub's  process  has  not  yet  been  made  avail¬ 
able  to  the  author. 

This  approach  has  been  successful  though  by  no  neans  perfected  or  completed.  Instead  of 
taking  a  mathematical  representation  of  the  velocity  structure,  a  method  was  devised  of  sampling 
the  structure  and  allowing  the  rays  to  traverse  the  samples  in  a  step-by-step  fashion.  As  com¬ 
puter  storage  is  often  limited,  a  method  was  devised  to  use  an  arbitrary  portion  of  the  data, 
save  the  information  of  the  emerging  rays  and  their  travel  times  and  amplitudes,  read  in  (or 
construct)  the  velocity  data  for  the  succeeding  portion,  and  continue  the  computations  through 
the  succeeding  portion.  The  number  of  portions  is  not  limited. 
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Multiple  reflections  were  included  by  continuing  with  a  reflected  ray  while  saving  the 
information  of  the  ray  that  was  incident  on  the  interface.  After  tracing  the  reflected  rays,  the 
program  returns  to  the  preceding  incident  rays  and  continues  tracing  through  the  interface.  As 
many  as  ten  multiple  reflections  have  been  allowed  for  in  the  exploratory  programs,  but,  again, 
this  number  can  be  as  high  as  wished,  and  can  be  terminated  either  by  limiting  the  number  or 
monitoring  the  amplitude  upon  reflection. 

An  extension  to  three-dimensional  ray  tracing  was  developed.  For  the  purpose  of  investi¬ 
gation,  a  simple  10  x  10  x  10  grid  was  employed  with  good  results. 

Thus  far  a  simple  approximation  has  been  employed.  Only  velocity  data  has  been  used. 

Snell's  law  is  employed  with  Fresnel  reflection  coefficients.  The  amplitude  has  been  computed 
by  fixed  exponential  decrement  and  spherical  spreading  where  the  initial  amplitude  is  divided 
by  the  distance  propagated. 

However,  the  use  of  these  simple  computations  are  considered  preliminary.  Extension  to 
more  accurate  seismological  approximations  will  naturally  follow.  With  inclusion  of  density 
and  S-wave  data,  Zoeppritz's  equations  can  be  implemented.  As  the  travel  times  are  computed, 
diffractive  effects  can  be  investigated  for  given  frequencies.  Ground  roll  can  be  added.  Syn¬ 
thetic  seismograms  may  be  possible.  Although  it  does  not  appear  that  propagation  computed 
with  the  wave  equations  will  naturally  follow  as  an  extension  of  this  method,  Keller's  geometric 
diffraction  theory  may.  Also,  by  dynamic  programming,  changing  of  the  velocity  structure 
with  iterations  appears  to  be  a  natural  extension  of  the  method. 

Previously,  an  IBM  7090  computer  and  Calcomp  plots  were  used  for  batch  processing.  The 

University's  time-sharing  IBM  360/67  computer  has  recently  become  operational,  and  a  change- 
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over  is  being  made  to  this  facility.  This  computer  has  about  2x10  bytes  of  storage.  Testing 
programs  have  been  run  on  the  360,  for  example,  for  development  of  the  equations  for  ray  tracing 
in  any  direction  (360°  freedom)  with  any  slope.  CRT  output  is  available  with  the  new  computer. 

The  batch  processing  on  the  7090  was  slow,  as  a  minimum  24-hr  turnaround  time  was  required. 
Difficulties  have  been  encountered  in  the  changeover,  and  during  some  periods  the  program  has 
seemed  to  retrogress,  which  is  all  too  common  with  computer  program  development.  When 
new,  more  demanding  data  were  added,  deficiencies  have  successively  been  discovered.  For 
example,  channelizing  through  successive  critical  reflections  due  to  the  sampling  interval  and 
computing  the  first  segment  after  reflection  have  been  recent  problems. 

The  author  presented  a  paper  on  this  method  in  April  1968  [13]. 

Preliminary  results  are  shown  in  figures  41-44.  These  are  not  presented  as  finished  prod¬ 
ucts  but  as  illustrations  of  the  potential  of  the  method.  Program  glitches  have  been  deliberately 
left  in  to  give  an  idea  of  some  of  the  problems  which  have  been  encountered.  The  average  execu¬ 
tion  time  during  computation  of  the  figures  has  been  30  sec,  with  plotting  times  of  about  2  min.  Of 
particular  significance  are  figure  45,  where  the  Jeffrey- Bullen  p-wave  travel- time  curve  is  closely 
duplicated  with  only  40  samples  of  the  radial  velocity  profile,  and  figure  43,  where  an  isometric 
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projection  of  the  results  of  three-dimensional  ray  tracing  is  shown.  For  the  case  of  the  spherical 
earth,  the  power  in  this  method  lies  in  the  capability  of  introducing  lateral  inhomogeneities. 

When  used  for  structural  determination  with  regional  disturbances,  this  ray-tracing  tech¬ 
nique  can  be  used  as  shown  in  figure  46.  For  local  structural  determination  from  teleseisms, 
it  can  be  used  as  shown  in  figure  47. 

7.2.  CURRENT  PROGRESS 

Currently,  FORTRAN  IVG  is  being  used  on  the  360  computer  for  two-dimensional  ray  trac¬ 
ing  with  complete  ray-direction  and  slope  freedom.  (In  MAD  language,  for  the  7090  computer, 
three-dimensional  ray  traces  were  made  with  a  representation  of  a  simplified  structure.)  A 
"patchwork"  program,  in  which  data  can  be  successively  read  in  after  tracing  through  and  re¬ 
taining  values  from  the  preceding  data,  has  been  perfected  and  the  resulting  travel-time  curves 
produced  by  the  emerging  rays  have  been  drawn  by  the  computer.  This  program  will  obviate  the 
need  for  large  memory  for  detailed  ray  tracing.  A  program  for  Including  head  waves  has  been 
written.  A  method  of  using  Zoeppritz'  equations  has  been  checked  out  and  is  the  next  planned 
addition  to  the  basic  program. 

7.3.  DISCUSSION 

The  program  in  the  appendix  will  reproduce  most  of  the  waves  produced  in  the  half-earth 
(with  the  exception  of  the  diffracted  P),  such  as  P,  PP,  PcP,  PKP  in  all  branches,  PKKP,  PKiKP. 
With  minor  modification  PKKP  can  be  obtained.  When  this  program  is  tested  with  102  P-wave 
velocity  values  sampled  from  the  Jeffrey- Bullen  velocity  vs.  depth  graph,  the  travel  times  for 
the  various  waves  are  usually  within  1%  of  those  found  in  the  Jeffrey- Bullen  travel-time  tables. 
The  errors  are  due  to  the  inaccuracies  of  taking  velocity  values  from  a  small  graph,  the  margin 
of  error  in  the  original  formulation  of  the  velocity  values,  and  the  result  of  sampling  a  relatively 
few  values.  In  ray  tracing,  any  error  due  to  sampling  in  the  initial  portion  of  the  rays  tends  to 
become  exaggerated  by  a  'lever-arm"  effect.  The  length  of  the  two  legs  of  PP  often  differ  by 
several  percent,  but  the  travel  times  are  more  accurate  due  to  the  fact  that  they  are  propagating 
through  the  proper  velocity  structure.  The  available  memory  of  the  computer,  of  course,  dictates 
the  amount  of  data  which  can  be  handled  at  one  time.  For  the  program  in  the  appendix,  all  the 
data  are  accepted  at  one  pass. 

In  the  computer  runs  with  102  samples  to  the  earth's  center,  the  average  time  to  trace  each 
ray  and  obtain  the  travel  time  was  1  sec  of  central  processing  unit  time. 

At  the  present  time  graphical  output  has  not  yet  been  incorporated  into  the  computer  using 
the  program  listed  in  the  appendix. 


FIGURE  41.  ILLUSTRATION  OF  SEISMIC  DATA  DISPLAY.  20  x  80  sample  points  read  in 

20  x  40  "patches." 
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TRAVEL  TIMES  DISTANCES  AMPLITUDES 
2.83  4.49  9.24 
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FRONT  VIEW  (I,J) 


SIDE  VIEW  (I,K) 


THREE-DIMENSIONAL  IOXIOXIO  SAMPLES 


FIGURE  43.  THREE-DIMENSIONAL  RAY  TRACING  THROUGH  A  VERY  SIMPLE  CUBE  WITH 

HORIZONTAL  LAYERS 


ARRAY  AREA 

SHOT  ^  '' 

POINT  SURFACE 


FIGURE  46.  RAY  TRACING  FOR  REGIONAL  SHOTS.  Travel  times,  amplitudes,  and 

angle  of  emergence  also  computed. 
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ARRAY  AREA 


INITIALIZATION  OF  RAYS  FROM  EACH  SAMPLING  POINT 
INITIAL  RAY  DIRECTIONS  FROM  SPHERICALLY  SPREADING  WAVES 


FIGURE  47.  RAY  TRACING  FOR  TELESEISMS.  Travel  times,  amplitudes,  and  angles 

of  emergence  also  computed. 
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Appendix 

SEISMIC- RAY  TRACING  BY  DIGITAL  MODELING 
P.  L.  Jackson 

Seismic  rays  can  be  traced  through  velocity  and  density  structures  of  any  degree  of  in¬ 
homogeneity.  A  digital  computer  is  used  to  trace  and  graphically  present  the  rays,  and  to  com¬ 
pute  the  travel  times  and  amplitudes  of  those  emerging  at  the  surface.  The  number  of  multiple 
reflections  can  be  selected  or  can  be  limited  by  a  minimum  amplitude  upon  reflection.  A  sam¬ 
pling  of  the  structure  on  two-  or  three-dimensional  grids  comprises  the  data  for  the  ray  tracing. 
As  the  data  can  be  "patched"  together  and  the  continuation  of  rays  traced  through  each  data 
"patch,"  only  a  small  amount  of  computer  memory  is  required;  thus  detail  is  not  limited  by  com¬ 
puter  memory.  Calcomp  plots  of  ray  traces  through  several  representative  velocity  structures 
have  been  made.  Included  is  an  earth  quadrant  in  which  the  Jeffrey- Bullen  P-wave  velocity 
profile  was  used  for  data.  The  program  can  be  extended  to  include  more  techniques  of  seismic 
interest,  such  as  Zoeppretz's  equations,  head  waves,  Rayleigh  waves.  One  goal  of  this  work  is 
to  present  the  output  on  a  CRT  so  that  the  results  of  velocity  and  density  data  changes  can  be 
immediately  perceived. 

Because  ray  tracing  for  a  spherical  earth  is  the  only  present  program  written  in  FORTRAN 
IVG,  only  this  program  will  be  described  in  this  report.  This  program  displays  the  basic 
approach  of  the  ray- tracing  procedure,  and  can  be  expanded  or  altered  by  the  reader.  The 
program  and  the  flow  diagram  (fig.  48)  follow  the  appendix. 

Snell's  law  is  used  to  iteratively  penetrate  through  an  array  (in  this  description  a  two- 
dimensional  array).  Each  element  of  the  array  is  a  velocity  value  representing  the  velocity  at 
a  depth  corresponding  to  the  row  of  the  array  and  a  distance  corresponding  to  the  column  of  the 
array.  A  corresponding  array  represents  the  slope  found  at  each  velocity  value.  As  the  ray 
penetrates  a  fixed  distance  at  an  angle  determined  by  Snell's  law,  the  velocity  for  the  incident 
value  and  slope  value  nearest  to  the  extent  of  penetration  are  taken  to  compute  the  angle  of  the 
next  penetration  of  the  fixed  distance  through  the  array.  A  "probing"  extension  of  the  fixed 
distance  at  the  preceding  angle  is  made  to  determine  the  velocity  of  the  emerging  region.  The 
velocity  of  the  incident  region,  the  angle  of  incidence,  and  the  velocity  of  the  emergent  region 
are  thus  determined  for  a  succeeding  application  of  Snell's  law.  If  desired,  a  change  of  one  set 
of  indices  can  be  made  to  use  the  slope  of  the  emergent  region  instead  of  the  incident  region. 

As  each  fixed  distance  of  the  penetrating  ray  is  incremented,  this  distance  is  divided  by  the 
velocity  value  through  which  the  ray  is  penetrating  to  compute  the  travel  time  for  this  incre¬ 
ment.  These  segmental  travel  times  are  successively  summed  so  that  upon  emergence  at  the 
surface  the  total  travel  time  is  retained. 
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When  the  ray  "probes"  into  an  emerging  region  whose  velocity  differs  by  a  specified  amount 
and  in  a  region  which  can  be  limited,  reflection  occurs.  The  value  of  the  reflecting  angle  is 
determined  by  the  angle  of  incidence,  and  the  penetration  continues  with  the  initialized  angle  of 
reflection.  In  case  of  critical  reflection  no  values  are  retained  from  the  reflecting  interface. 

With  noncrltical  reflection,  the  value  of  the  incident  angle  is  retained,  along  with  the  summed 
travel  times  to  the  interface.  After  the  reflected  ray  emerges  or  exceeds  the  bounds  of  the 
array,  control  is  returned  to  the  incident  ray  so  that  it  can  penetrate  beyond  the  interface.  The 
number  of  multiple  reflections  can  be  specified  or,  if  desired,  terminated  if  the  amplitude  is  suf¬ 
ficiently  small,  although  this  latter  condition  is  not  in  the  accompanying  program. 

When  any  ray  is  sufficiently  close  to  the  representation  of  the  surface  of  the  earth,  its 
azimuthal  distance  and  its  length  (in  kilometers),  the  travel  time,  and  the  amplitude  are  printed 
out  before  a  new  direction  from  the  source  or  returning  to  a  reflection  point  is  initiated  for 
continuation  of  the  incident  ray. 

The  program  contains  two  subroutines.  The  first,  entitled  "WORLD,"  accepts  radial  velocity 
samples  to  any  depth  and  constructs  the  tw>  dimensional  velocity  and  slope  arrays  to  correspond 
to  a  curved  earth.  When  the  half- earth  is  used,  it  constructs  the  array  for  half  of  the  earth.  T'u 
second,  entitled  "SINDET,"  accepts  the  angle  from  the  vertical,  the  velocity  value  of  the  incident 
region  and  its  slope,  and  the  velocity  value  from  the  emerging  region.  It  then  computes  the  emerg¬ 
ing  angle  by  means  of  Snell's  law.  Any  angle  from  the  vertical  from  0°  to  360°  and  any  slope  from 
-90°  to  +90°  is  accepted  by  SINDET.  In  case  of  angles  which  are  critical  or  beyond,  this  informa¬ 
tion  is  returned  to  the  main  program  so  that  a  reflection  will  take  place. 

Because  a  rectangular  array  is  used,  the  conversion  to  flat  earth  is  straightforward.  The 
subroutine  "WORLD"  is  not  used  for  flat  earth,  and  the  rays  emerge  at  the  first  row  of  the 
array.  Actually  the  flat  earth  is  more  economical  to  compute  because  an  arctangent  computation 
is  not  required  for  each  element  of  the  array. 
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MAIN  PROGRAM  SUBROUTINE  WORLD  AND  SUBROUTINE  SINDET 


> 

1 

DIMENSION  A( 203* 102),BC203, 102 ) , PI C 10 >, P JC 10 > , 

> 

2 

1CI J( 10),0N( 10) « ANT ( 1 0) »D( 102)»  AT  (  10) 

> 

3 

ICRIT*0 

> 

4 

NAMELI ST/NAM  1 /NO IN>  DELN»BEGIN»NMULT»ORI GV«0RI GH«H ‘ 

> 

5 

1, I SAMP* IPORSH*  JPORSH*  DIFRAC 

> 

6 

READC4.NAM 1 ) 

> 

7 

REA DC  5*  500 ) C  DC  I > , I -  1 , IPORSH) 

> 

8 

CALL  WORLDC  D>  A>  B»  I  SAMP#  IPORSH*  JPOHSH) 

> 

9 

WRITEC 1, 600) ( (AC  I,  J),  I«l,203),  J-l,  102) 

> 

10 

WRITEC  1, 600 )C  CBC I»  J)>  1-1,203)* J- 1> 102) 

> 

1  1 

DATA  PIT/3.141  5927/ 

> 

18 

DATA  DEC/-. 00003/ 

> 

13 

C 

> 

14 

C 

17  CALCOMP  CARDS 

> 

15 

c 

> 

16 

DO  20  L- 1 , NOIN 

> 

17 

IFCBEGIN.GT. 1.57)  CALL  SYSTEM 

> 

18 

III-O 

> 

19 

WRI TEC  7, 50 1 )  BEGIN 

> 

20 

ANGLE2-BEGIN 

> 

81 

PNM-1000. 

> 

28 

DO  30  1-1,  10 

> 

23 

ATCI >=0. 

> 

24 

30 

ONC I )-0. 

> 

25 

IVW-0 

> 

26 

IZ-1 

> 

27 

ANTC IZ)- 1 . 

> 

28 

PICIZ)-ORIGV 

> 

29 

P JC I Z)-ORI GH 

> 

30 

I-PI C I Z ) ♦• 500  1 

> 

31 

J-PJC I Z  > ♦. 500 1 

> 

32 

40 

IFCIZ)  20,20,60 

> 

33 

60 

IU-I 

> 

34 

IV-J 

> 

35 

ANGLE1-ANGLE2 

> 

36 

C 

> 

37 

C 

2  CALCOMP  CARDS 

> 

38 

C 

> 

39 

I F( ANGLE1 . GT. 7 .  )  GOTO  52 

> 

40 

G-H*ACI, J)/13.62 

> 

41 

HHH«G*COS( ANGLE1 ) 

> 

42 

GGG- G*S INC  ANGLE 1 ) 

> 

43 

PI(IZ)-PICIZ) +HHH 

> 

44 

PJC  IZ)-P JC IZ)  +  GGG 

> 

45 

ATCIZ)-ATCIZ)+DI FRAC  *G/(A(I,J)*60.) 

> 

46 

ONC IZl-ONC IZ)+DIFRAC*G 

> 

47 

PI  I- 1 SAMP- PI ( I Z) 

> 

48 

PJJ-PJCIZ)-1. 

> 

49 

RAD-SQRTCPI I*PI I+PJJ*PJJ> 

> 

50 

IFCRAD- I SAMP+1 . >  70,70,80 

> 

51 

80 

IFCBCI, Jl.LE.O. )  DISTA— 57.2957795*BC  I,  J) 

> 

52 

IFCBCI, Jl.GT.O.  )  DI STA-57.2958*C PI T-BC I , J>  > 

> 

53 

AMPL-  CPNM*ANT(  IZ)/ON(IZ)  >*EXPC  DEC*ON(  I Z)  ) 

> 

54 

DKM-11 1 . 4726*DI STA 

> 

55 

WRITEC 7 ,501)  DISTA, DKM,AT(IZ), AMPL 

> 

56 

III-III+l 

> 

57 

IFCIZ. EQ.l. AND. ICRIT.EQ.O)  GOTO  210 

> 

58 

GOTO  52 

> 

59 

70 

IF(PI(IZ).GT.. 5. AND.  P  J(  I Z)  • GT. • 5 .AND. PI (IZ).LT. 202.5 

> 

60 

1  .AND.PJC IZ) .LT. 101.5)  GOTO  90 

> 

61 

52 

IVW-0 

> 

62 

ANTC  IZ-1  >-C 1-ANTC I Z) )*ANTC I Z- 1 ) 

> 

63 

IZ-IZ-1 

> 

64 

IFCIZ)  20,20,100 
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MAIN  PROGRAM  SUBROUTINE  WORLD  AND  SUBROUTINE  SINDET  (Continued) 


> 


> 

> 

> 

> 

> 


65 

100 

I  =  P I t I Z)  +  <  5 

66 

J-PJt IZ1+. 5 

67 

ANGLE.2-C1  JtIZ) 

68 

GOTO  40 

69 

90 

I-PItIZ)*.5 

70 

J*PJt IZ)*. 5 

71 

C 

3  CALCOMP  CARDS 

72 

C 

73 

PlI*PItIZ)+HHH 

74 

PJJ»PJ( 1Z1+GGG 

75 

IM«PI I+.5 

76 

IS-PJJ+. 5 

77 

PI I-I SAMP-PI  I 

78 

PJJ*PJJ-1. 

79 

C********C0MPUTE  RADIUS******** 

80 

RAD»SQRT  <  P I I*PII+PJJ*PJJ> 

81 

WRITE! 1, 5011  PI ( I Z), PJt IZ),AtI,J),BtI,J) 

82 

C* ****** *CHECK  FOR  BEING  IN  BOUNDARIES******** 

83 

IFtRAD.GT. ISAMP- 1 )  GOTO  80 

84 

I  Ft IM.GT.2*IP0RSH-  1  .OR. IM.LT. 1 .OR. IS. GT. JPORSH 

85 

1. OR.IS.LT. 1)  GOTO  52 

86 

IF<A< IM, IS)-Af I. J) )  120,40.120 

87 

C********CALL  SUBROUTINE  TO  COMPUTE  DIRECTIONS******** 

88 

120 

CALL  SINDET<ANGLE1,ANGLE2,AC I, J) , At IM, I S ) ,  Bt  I  » 

89 

1 J), ICRIT) 

90 

I  Ft ICRIT+IUW.GE.2)  GOTO  40 

91 

IFtICRIT.GE.l)  GOTO  230 

92 

130 

IF(ABS(A< I, J)-At IM. IS) ) .GE. 1 . 2 . AND. I Z. LE. NMULT 

93 

l.AND.IUW.LT.l .AND. At IM, IS). GT. 3.. AND.HAD.LT. 64. 

94 

2)  GOTO  135 

95 

IUW-0 

96 

GOTO  40 

97 

135 

IZ-IZ+1 

98 

0N<IZ)-0N<IZ-1 ) 

99 

AT(IZ)-ATUZ-l) 

100 

PItIZ)*PIt IZ-1 ) 

101 

PJt  I  Z)"PvJ<  I  Z-  1  ) 

102 

ANT< IZ)-ANT( IZ- 1 )*(AC I, J)-A( I U, IU) )/CA( I , J)*A( I U, IV) ) 

103 

Cl J< IZ-1 )“ANGLE2 

104 

230 

I  Ft  I Z-NMULT)  210,220,220 

105 

220 

IZ-IZ-1 

106 

GOTO  40 

107 

210 

ANGLE2-PIT-ANGLE1*2.*B( I ,  J) 

108 

I Ft ANGLE2.GE.2.*PIT)  ANGLE2- ANGLE2-2 . *PI  T 

109 

IFtANGLE2.LT.O.)  ANGLE2- ANGLE2*2. *PI T 

1  10 

IFtIII.GT.1)  GOTO  52 

1  1  1 

IVW-1 

112 

GOTO  40 

113 

20 

BEGIN-BEGIN*DELN 

114 

C 

115 

C 

1  CALCOMP  CARD 

116 

C 

117 

500 

FORMATt 10F4.2) 

118 

501 

FORMAT t  6F10 . 2 ) 

119 

700 

FORMAT  t  F3. 2 ) 

120 

600 

FORMATt '  * 10F5 .2  ) 

121 

CALL  SYSTEM 

122 

END 

123 

C 

124 

C 

125 

SUBROUTINE  WORLD t  D, VEL, SLP, I  SAMP, I PORSH,  JPORSH) 

126 

DIMENSION  Dt 1 02 ) , VELf 203, 1 02) , SLPt 203,  102) 

127 

PPI-3. 1415926/2. 

128 

DO  10  1*1, IPORSH 

129 

SLPt IPORSH, I )*PPI 

130 

SLPt I , l)-0. 

MAIN  PROGRAM  SUBROUTINE  WORLD  AND  SUBROUTINE  SINDET  (Concluded) 


> 

131 

VELI IPORSH# I >«  DC IP0RSH+1-I) 

> 

132 

10 

VELI I* 1 >-D(  I  > 

> 

133 

IPOR- IPORSH- 1 

> 

134 

DO  20  I«  1 » IFOR 

> 

135 

DO  20  J-2#  J.-’ORSH 

> 

136 

P- IPORSH- I 

> 

137 

Q«  J-l 

> 

138 

POP*P+Q*Q 

> 

139 

H-SQRTIPQ) 

> 

140 

IFCH.GT. IP0R+.5001)  GOTO  20 

> 

141 

R-P/Q 

> 

142 

THETA-ATANCR) 

> 

143 

SLPI I, J5-THETA-PPI 

> 

144 

K«  IPORSH-H+. 5001 

> 

145 

VELI I  #  J) ■ Dt  K) 

> 

146 

20 

CONTINUE 

> 

147 

IF(IPORSH.NE.ISAMP)  GOTO  40 

> 

148 

DO  30  1“  1 »  IPOR 

> 

149 

DO  30  J=  1 »  JPORSH 

> 

150 

VELI I  +  IPORSHi J)  =  VEL( IPORSH-I,  J) 

> 

151 

30 

SLPI I  +  IPORSHi J)«-SLP( IPORSH- 1 » J) 

> 

152 

40 

RETURN 

> 

153 

END 

> 

154 

C 

> 

155 

C 

> 

156 

SUBROUTINE  S INDETt Cl U#CI  A I J, AMS# BI I CRI T ) 

> 

157 

ICRIT-0 

> 

158 

DATA  PI/3.1415927/ 

> 

159 

DATA  PPI/1 .57579631/ 

> 

160 

AA= AMS/A I J 

> 

161 

GH=0  • 

> 

162 

HG«1. 

> 

163 

CC-CIU-BI J 

> 

164 

IF(CC.EQ.0..0R.CC.EQ.PI.0R.CC.EQ.2.*PI )  GOTO 

> 

165 

IFICIU.GT.O. . AND. Cl U.LE.PI )  GOTO  101 

> 

166 

IFICC.LT. 3. *PP I )  GOTO  112 

> 

167 

GH-2.*PI 

> 

168 

GOTO  1000 

> 

169 

101 

IFICC.GT.PPI )  GOTO  112 

> 

170 

GOTO  1000 

> 

171 

112 

IH'PI 

> 

172 

HG--1. 

> 

173 

1000 

FFF« AA*S INI CC ) 

> 

174 

IFIABSIFFF) . GE. 1 . )  GOTO  13 

> 

175 

ARC-HG*ARSINIFFF) 

> 

176 

Cl  J-ARC+BI J*GH 

> 

177 

IFICI J.GT.2.+PI )  Cl J»CI J-2.*PI 

> 

178 

IFICI J.LT.O. )  Cl J»CI J*2.*PI 

> 

179 

GOTO  14 

> 

180 

13 

ICRIT-1 

> 

181 

GOTO  14 

> 

182 

1 1 

Cl J»CIU 

> 

183 

14 

RETURN 

> 

184 

END 

# 

END  OF 

FILE 

# 
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DIMENSION  A,B,P1, 

P.1.CIJ,TT,AT)ON) 

ANT,D 


DATA  INPUTS 
ICR1T-0,  PIT-TT 
DEC. 

NAMELIST  DICTIONARY 
NOIN.DELN,  BEGIN, 

NMU  LT,IPORSH,SAMP, 
OR1GV.ORIGH.H, 

DIFRAC 


READ  VELOCITY- 
DEPTH  SAMPLES:  D 


X 


CALL  WORLD  (D, 
A,B,SAMP,lPORSH,JPORSH) 

I 


WRITE  VELOCITY 

AND  SLOPES:  A,B 

L- NOIN  - - 

- - 

L-l  " 

do  o 

3 

ET 

BEGIN>  1.57 


NO 


WRITE  INITIAL 
ANGLE:  BEGIN 


ANGLE2-  BEGIN 
1  .  " 


PNM-1000 


J 


J 

T 

o 

n 

i-i 

AT(I)-0 

Noa)-o 

IVW-0 

IZ-I 

ANT(IZ)- 1 
PI(IZ)-ORI  GV 
PJ(IZ)-ORIGH 
I-PI(IZ)+.5001 
J-PJ(IZ)+.500I 


NO 

IU-I 
IV-J 
ANGLE  1- 

-ANGLE2 

ANGLE  1>  7.0 


YES 


G*HxA(I,J)/13.62 
HHH-GxCOS  (ANGLE  1) 
GGG-GxSIN  (ANGLE1) 
PI(IZ)-PI(IZ)+HHH 
PJ(IZ)-PJ(IZ)+GGG 
AT(IZ)-AT(IZ)+DIFRAC 
xG/A(I,J) 

ON  (IZ )-  ON  dZ  )+DI  FRAC 
xG 

PU-SAMP-PI(IZ) 

PJJ-PJ(IZ)-I. 

RAD-/PH2+PJJ2 


FIGURE  48.  FLOW  CHART  FOR  MAIN  PROGRAM  SUBROUTINE  WORLD  AND  SUB¬ 
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FIGURE  48.  FLOW  CHART  FOR  MAIN  PROGRAM  SUBROUTINE  WORLD  AND  SUB¬ 
ROUTINE  SINDET  (Continued) 
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MAIN  (PAGE  3) 


CALL  SINDET 
(ANGLE  1,  ANGLE2, 
A(I,J),  AdM.IS), 
Ba,J),ICFUT 


1Z-1Z-0 

ONflZ)-ONdZ-l) 

ATdZJ-ATdZ-l) 

piazj-  piaz-i) 

pjaz)-pj(iz-u 

ANT(IZ)-ANTOZ-I) 

x(AflM,IS)-Aa,J))/ 

(AaM,is)+Aa,j)» 

CU(IZ-1)-  ANGLE2 


FIGURE  48.  FLOW  CHART  FOR  MAIN  PROGRAM  SUBROUTINE  WORLD  AND  SUB¬ 
ROUTINE  SINDET  (Continued) 
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FIGURE  48.  FLOW  CHART  FOR  MAIN  PROGRAM  SUBROUTINE  WORLD  AND  SUB¬ 
ROUTINE  SINDET  (Continued) 
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START 


FIGURE  48.  FLOW  CHART  FOR  MAIN  PROGRAM  SUBROUTINE  WORLD  AND 

ROUTINE  SINDET  (Continued) 
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